/*
 *  polynomial.cpp
 *  radici
 *
 *  Created by Nicola Gigante on 04/03/08.
 *  Copyright 2008 Nicola Gigante.
 *
 */

#include "polynomial.h"

const int rootsMinIterations = 1;
const Real dEpsilon = 1E-18;

Polynomial::Polynomial(const Polynomial &p)
{
	m_terms = p.m_terms;
	m_complex = p.m_complex;
	m_variable = p.m_variable;
}

Polynomial::Polynomial(Complex known_term)
{
	m_terms.append(known_term);
	m_complex = false;
	m_variable = 's';
}

Polynomial::Polynomial(Real known_term)
{
	m_terms.append(known_term);
	m_complex = false;
	m_variable = 's';
}

/**
 Questa funzione statica crea un polinomio partendo da un numero complesso.
 Il risultato è diverso dal costruttore analogo.
 Infatti mentre il costruttore Polynomial(Complex) imposta il numero complesso come
 termine noto, questo metodo crea un polinomio P(j) = a + j b partendo dalle componenti
 cartesiane del polinomio stesso.
 \see Polynomial(Complex)
 */
Polynomial Polynomial::complexToPolynomial(Complex c)
{
	Polynomial poly;
	poly.addTerm(0, real(c));
	poly.addTerm(1, imag(c));
	poly.setVariable('j');
	poly.m_complex = true;
	
	return poly;
}

void Polynomial::setTerm(int exp, Complex k)
{
	if(exp > grade())
	{
		// Fill empty terms
		for(int i = grade(); i < exp; i++)
			m_terms.append(0);
	}
	m_terms[exp] = k;
}

void Polynomial::addTerm(int exp, Complex k)
{
	if(exp > grade())
	{
		// Fill empty terms
		for(int i = grade(); i < exp; i++)
			m_terms.append(0);
	}
	m_terms[exp] = m_terms.at(exp) + k;
}

Complex Polynomial::term(int exp) const
{
	if(exp > grade())
		return 0;
	else
		return m_terms.at(exp);
}

long double Polynomial::maxTermModule()
{
	long double m = 0;
	
	foreach(Complex term, m_terms)
		if(abs(term) > m)
			m = abs(term);
	
	return m;
}

/**
 Questo metodo è utilizzato principalmente all'interno dell'algoritmo di ricerca
 delle radici. Ogni coefficiente del polinomio viene diviso con il coefficiente
 del termine di grado massimo. Il risultato è un polinomio con il termine di grado
 massimo uguale a 1. Le radici del polinomio non vengono alterate
 */
Polynomial Polynomial::norm()
{
	Polynomial poly;
	
	Complex c = term(grade());
	
	for(int i = 0; i < grade(); i++)
	{
		poly.setTerm(i, term(i)/c);
	}
	
	poly.setTerm(grade(), 1);
	
	return poly;
}

Complex Polynomial::eval(Complex val)
{
	Complex r = 0;
	
	for(int i = 0; i <= grade(); i++)
		r = r + m_terms[i]*pow(val, i);
	
	return r;
}

/**
 L'addizione viene effettuata addizionando tra loro i coefficienti
 di pari grado, come dalle regole dell'algebra.
 \note La funzione accetta un parametro di tipo Polynomial ma grazie al costruttore
	   di conversione Polynomial(Complex) e Polynomial(Real) è possibile addizionare
	   anche solo un numero complesso o reale, di fatto variando il termine noto.
 */
Polynomial Polynomial::add(Polynomial p)
{
	Polynomial poly = *this;
	
	for(int i = 0; i <= p.grade(); i++)
		poly.addTerm(i, p.term(i));
	
	return poly;
}

/**
 La moltiplicazione viene effettuata come dalle regole dell'algebra, moltiplicando
 ogni termine del primo polinomio per ognuno dei termini dell'altro polinomio.
 */
Polynomial Polynomial::multiply(Polynomial p)
{
	Polynomial poly;
	
	for(int i = 0; i <= this->grade(); i++)
		for(int j = 0; j <= p.grade(); j++)
			poly.addTerm(i + j, this->term(i) * p.term(j));
	
	return poly;
}

/**
 Moltiplicazione con un numero complesso.
 La moltiplicazione viene effettuata come dalle regole dell'algebra, moltiplicando
 ogni termine del polinomio per il numero complesso c.
 */
Polynomial Polynomial::multiply(Complex c)
{
	Polynomial poly;
	
	for(int i = 0; i <= this->grade(); i++)
		poly.addTerm(i, this->term(i) * c);
	
	return poly;
}

/**
 Moltiplicazione con un numero reale.
 La moltiplicazione viene effettuata come dalle regole dell'algebra, moltiplicando
 ogni termine del polinomio per il numero reale r.
 */
Polynomial Polynomial::multiply(Real r)
{
	return this->multiply(Complex(r));
}

Polynomial Polynomial::reverse()
{
	Polynomial p = *this;
	
	::reverse(p.m_terms.begin(), p.m_terms.end());
	
	return p;
}

/**
 Questo metodo è di gran lunga il più importante e utilizzato della classe.
 Ricerca le radici complesse del polinomio tramite il metodo di Durand-Kerner.
 Durand-Kerner è un metodo ad approssimazioni successive.
 Si è riscontrato nella pratica che in presenza di errori di approssimazione
 nelle operazioni sui numeri in virgola mobile il metodo potrebbe non convergere
 causando un loop infinito. Questa situazione, comunque rarissima, viene intercettata
 e il metodo lancia un eccezione di tipo QString con valore "Errore di convergenza. Reset".
 La spiegazione dettagliata del metodo di Durand-Kerner è disponibile nella documentazione
 tecnica del programma.
 */
QList<Complex > Polynomial::roots()
{	
	QList<Complex > roots;
	QList<Complex > prev;
	
	// Max coefficient must be 1
	Polynomial poly = norm();
	
	// Compute initial guessed roots
	srand(time(NULL));
	Real m = poly.maxTermModule();
	for(int i = 0; i < grade(); i++)
	{
		Complex r = polar((double)(rand() * m / RAND_MAX), rand() * 2*M_PI / RAND_MAX);
		if(imag(r) != 0 && abs(r) != 1)
			prev.append(r);
		else
			i--;
	}
	
	roots = prev;
	
	bool changed = false;
	Complex num, den;
	bool good = false;
	int iters = 0;
	do
	{
		changed = false;
		
		if(++iters == 1000)
			throw QString("Errore di convergenza. Reset");
		
		for(int i = 0; i < prev.count(); i++)
		{
			num = poly.eval(prev[i]);
			den = 1;
			for(int j = 0; j < prev.count(); j++)
			{
				if(i != j)
					den = den*(prev[i] - prev[j]);
			}
			
			roots[i] = prev[i] - num/den;
			
			//if(roots[i] != prev[i])
			if(!complexEqual(roots[i], prev[i]))
				changed = true;
		}
		if(!changed)
			good = true;
		
		prev = roots;
	}while(!good);
	
	return roots;
}


void Polynomial::print()
{
	printf("g: %d\n", grade());
	for(int i = grade(); i >= 0; i--)
	{
		printf("%f%s%s", m_terms[i].real(), i > 0 ? "x" : "", i > 1 ? "^" : "");
		if(i > 1)
			printf("%d", i);
		if(i > 0)
			printf(" + ");
	}
	printf("\n");
}

